arXiv:1503.07297vl [cs.CV] 25 Mar 2015 



ELSEVIER 


Available online at www.sciencedirect.com 

ScienceDirect 

Procedia Computer Science 00 (2015) 1-[40| 


Procedia 

Computer 

Science 


A Brief Survey of Recent Edge-Preserving Smoothing Algorithms 

on Digital Images 


Chandrajit Pal^ Amlan Chakrabarti^, and Ranjan Ghosh^ 

Corresponding author mail:palchandrajit@gmailcom 1. 

A.K.Choudhury School of Information Technology 1,2, Institute of Radio Physics and ElectronicsS, University of Calcutta, 92 A. P. C. Road, 

Kolkata:700 009, India. 


Abstract 

Edge preserving filters preserve the edges and its information while blurring an image. In other words they are used to smooth 
an image, while reducing the edge blurring effects across the edge like halos, phantom etc. They are nonlinear in nature. Exam¬ 
ples are bilateral filter, anisotropic diffusion filter, guided filter, trilateral filter etc. Hence these family of filters are very useful in 
reducing the noise in an image making it very demanding in computer vision and computational photography applications like de- 
noising, video abstraction, demosaicing, optical-flow estimation, stereo matching, tone mapping, style transfer, relighting etc. This 
paper provides a concrete introduction to edge preserving filters starting from the heat diffusion equation in olden to recent eras, an 
overview of its numerous applications, as well as mathematical analysis, various efficient and optimized ways of implementation 
and their interrelationships, keeping focus on preserving the boundaries, spikes and canyons in presence of noise. Furthermore it 
provides a realistic notion for efficient implementation with a research scope for hardware realization for further acceleration. 

Keywords’.Am^oiropic diffusion. Field Programmable Gate Array, bilateral filter, non local means, system generator. Very 
High Speed Integrated Circuit Hardware Description Language, FPGA-in-the-loop, Peak Signal to Noise Ratio, Partial Differential 
Equation, Look Up Table, Digital Signal Processor. 


1. Introduction 

From the very beginning there a rapid increasing interest in models and methods for discontinuous phenomena, 
both in the computer graphics and vision community. The main focus lies in the identification of discontinuities in 
data perturbed by noise. To be more specific, in image data where noise needs to be removed from image while 
preserving basic significant features like gradients, jumps, spikes, edges and boundaries. 

Methods dealing with region boundary preservation have a long tradition in imaging and in their state of self 
applied art. Each one have their promising contributions from various branches of mathematical perspective in their 
respective application domain. 

An evolutionary review study has been performed with detail analysis of every method starting from diffusion 
process O , nonlinear Gaussian filters lO, bilateral filter|[5l|7l[8l[TTl to its extension the trilateral filter lfT3ll4^ . 

Edge preserving filtering is a technique to smooth images while preserving edges. It all started in 1985 with the 
work of L.P Yaroslavsky (Tl, with his neighborhood filter which averages pixels with a similar grey level value and 
belong- ing to the spatial neighborhood. In 1990 Perona and Malik O worked on realizing semantically meaningful 
edges using a diffusion process. To overcome some of its limitation Aurich and Weule in the year 1995 made a 
breakthrough invention on edge preserving method using non-linear Gaussian filters. Their concept laid the foundation 
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Stone of the origin of bilateral filters. It was later rediscovered by Smith and Brady IH , and Tomasi and Manduchi m 
in the year 1998 gave it its current name ’the Bilateral filter’. Ever since it has found a widespread use in various image 
processing applications like video abstraction, demosaicing, optical-flow estimation, stereo matching, tone mapping, 
stylization, relighting and denoising. Its various advantages lies firstly with its non-iterative and simple formulation 
where each pixel is replaced by a weighted average of its neighborhood pixels making it adaptive to application 
specific ambiance. Secondly it depends on few parameters (spatial domain and intensity range kernels). Finally with 
the superb contributions of various efficient numerical schemes the filter can be computed very fast reasonably on 
large images also in real time environment with dedicated graphics hardware. 

A great deal of implemented papers E III 11113 uni El [la have been studied, which explores the various feature 
of the edge preserving filters (mainly bilateral filter) constituting its behaviour and gives a detail understanding of 
its strengths and weaknesses in various application environment domains, the various optimization techniques for 
further betterment which consequently lead to its extension the next generation trilateral filter |[T3l|43l. It provides 
stronger noise reduction and better outlier rejection in high-gradient regions, and it mimics the edge-limited smoothing 
behavior of shock-forming PDFs by region finding with a fast min-max stack. This study is supposed to provide a 
good insight to the designers among the selection of the design methods of their choice deciding upon the various 
numerical analysis compatible with their applications in hand. 

This paper is organized as follows. Section 2 presents the diffusion technique employed to detect the edge, namely 
the anisotropic diffusion, next realizing the same through nonlinear Gaussian filters, linear Gaussian filtering and the 
nonlinear extension to the bilateral filter and determining the computational cost, number of iterations followed by 
the various effects on it by changing its parameters. Section 3 discusses the various ways to efficiently implement 
the bilateral filter together with the variations of each method to achieve the other and also the several finks mainly 
of the bilateral filter to other frameworks employing the different methods to interpret it. Section 4 describes various 
challenging applications of the edge-preserving filters. Section 5 describes the various interrelationships among the 
edge preserving filters with an in depth mathematical analysis of each transition. Section 6 exhibits the various 
modifications, extensions and substitutes of the bilateral filter to the next generation trilateral filters. 


2. Edge preserving through nonlinear diffusion to trilateral filtering 

To introduce the concept of edge preserving we begin with our discussion with the nonlinear diffusion technique 
by Perona and Malik El- They proposed a nonlinear diffusion method for avoiding the blurring and localization 
problems caused by linear diffusion filtering. It aims at reducing image noise without removing significant parts of 
the image content, typically fines, edges or other details that are important for the interpretation of the image. 


2.7. Motivation leading to realizing the anisotropic diffusion filter 
Suppose we have an edge/not edge estimation function E such that 


E{x,yj) = V/(x,y,0 (1) 

where 


dl dl 

VI(x,yj) = (-,-) 
ox oy 


( 2 ) 


If (x,y) is a part of an edge then naturally little smoothing should be applied unlike otherwise. Fortunately the 
gradient of the brightness function acts as a detector which tells whether (x,y) is a part of an edge or not denoted by 
V7(x,y, t) = (|^, ^). Fet Q c denote a subset of the plane and 7(.,t):Q ^ M be a family of gray scale images 
therefore anisotropic diffusion is defined as 


— = div(c(( X, y, t)AI)) = VcV I c(x, y, t)AI 
ot 


( 3 ) 


where A and V denotes the laplacian and gradient respectively, div(...) is the divergence operator and c(x,y, 7) is 
the diffusion coefficient, c(x,y, 7) = g(\\E(x,yj)\\) controls the rate of diffusion i.e blurring intensity according to 
\\E(x,y, Oil and is the function of the image gradient so as to preserve the edges in the image. When c(x,y, t) is large 
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Figure 1. In the first part(l), the gradient vector is (0, n) while (n>0). In the second part(2), the gradient vector is (-m , m) while (m >0 ). These 
points where the norm of the gradient is high, could be treated as edge points, end therefore be applied less blurring. 


(x,y) is not a part of an edge else otherwise. The edge estimation function is E(x,y,t)=V/(x,y, t) correspondingly the 
two functions for the diffusion coefficient, the corresponding two functions for the diffusion coefficients are 

c(\\E\\) = e-m\/kfand 

c(||£||) = 1/(1 + ||£||A)^ (4) 

where the constant k controls the sensitivity to the edges and is either chosen experimentally or as a function of the 
image noise. 


If gradient is (a,b), then a perpendicular vector to it, can be (a, -b). A dynamic kernel is contracted along the 
direction of the normal, ending in an elliptical kernel. 

We have successfully implemented this algorithm after identifying the parallelism inherent with it on reconfigurable 
hardware and have attained a prominent acceleration ca. We have presented the FPGA implementation of an edge 
preserving anisotropic diffusion filter for digital images. The designed architecture completely replaced the convolu¬ 
tion operation and realized the same using simple arithmetic subtraction of the neighboring intensities within a kernel, 
preceded by multiple operations in parallel within the kernel. Our proposed hardware design architecture completely 
substituted standard convolution operation |l4l, required for evaluating the intensity gradient within the mask. We 
used simple arithmetic subtraction to calculate the intensity gradients of the neighboring pixels, which saved 9 multi¬ 
plication and 8 addition operations per convolution respectively by computing only a single arithmetic operation for 
each gradient direction as per equation 5. Thus, reducing a huge amount of computational time and operations to 
a linear computational complexity. The hardware implementation of the proposed design shows better performance 
in terms of throughput ( increase of 10% ), total power ( reduction of 24% ) and a reduced resource utilization in 
comparison with the related research works In the later sections we will show the relationship of this algorithm with 
its counterparts. 

2.2. Edge preserving diffusion through Non-Linear Gaussian filters 

Aurich and Weule in the year 1995 m proposed a new diffusion method for edge preserving smoothing of images. 
Unlike the anisotropic diffusion methodology it is based on a modification of the way the solution of the heat conduc¬ 
tance equation is obtained by convolving the Gaussian kernel with the initial data. Their method employs a simple 
non-linear modifications of the Gaussian filters thereby overcoming the lengthy iteration steps and convergence prob¬ 
lems. A sequential chain of non-linear filters combining the significant as well as advantageous aspects of both the 
anisotropic diffusion and robust filtering has been applied overcoming their drawbacks. 
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Figure 2. 


The desired direction of the dynamic kernel is effected by the local gradients direction V/(x, y, t) and it is perpendicular to it. 



a) Original b) Anisotropic diffusion 


Figure 3. 


fig ’a’ denotes the original figure and ’b’ denotes the denoised image after applying anisotropic diffusion i.e Perona Malik filter. 


2.2.1. The non-linear Gaussian filtering approach 

A simple linear Gaussian filtering convolution is defined by Gf=g*/ where / is the signal and g is the Gaussian 
function. The image signals are defined on a discrete bounded set p of pixels and as a result region boundary effects 
needs to be taken care off. The linear Gaussian filter using the position-dependent renormalization of the weights is 
defined by 


Ga/(P) = / Z ^-411^ - pmp) (5) 

p qep 

= f(p) + - pWXfiq) - fip)) 

P qep 

with ga-yt)=expi-f llcrl) and Np =Y,qep ga-^M - P\\)- 

In order to preserve edges one should avoid averaging across edges. Thus if q and p are at different sides of an 
edgef(q) should be disregarded while calculating the average over a neigbourhood of p. Usually, however, the location 
of the edges is not yet known, even worse: Edge-preserving is done in order to detect edges more easily. Therefore 
the valueis not completey omitted in the averaging process, but its weight is lowered when it is likely that there 
is an edge between q and p. How to estimate this probability? The difference \f(q)-f(p)\ needs to be looked upon. The 
bigger it is, the more probable it seems that there is an edge between p and q. This has been realised by multiplying 
~ P\\) additional weight i/zifiq) - f(p)). As because the noise distribution in images in real situations 

is bell-shaped i.e Gaussian like so it is suggestive to choose if/ = go-^ i.e also as Gaussian, thereby generating the 
non-linear Gaussian filter given by 

= m + / E - p\\)g^ym - mum - mi (6) 

P qep 
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Figure 4. shows the effect of the parameter 77 with 1.0 and 1.3 at the left and right figures respectively. Courtesy of O. 


/ 

Gf 


-4 -3 



■> 


Figure 5. shows the effect of sharpening of a blurred edge. Courtesy of 0. 


where 

^P = 'Z - PiDgaAm - fip))- (7) 

qep 

It is well known that the Gaussian filter is the basic solution of the heat conductance equation, therefore the non¬ 
linear Gaussian filter can be utilized to calculate some form of anisotropic diffusion. With an aim to further decrease 
the noise an extra parameter 77 has been added to the nonlinear Gaussian filter equation [^resulting in 

Ga,,a.Ap) = fip) + ^ Z - p\\)g<TAPq) - fipMfiq) - fip))- (8) 

p qep 

Experimental observations have revealed that larger the value of 77 narrower is the noise distribution of the filters as 
shown in fig 4 below. 

The value of the parameters namely (t^ and cr^ plays an important role in determining the shape of the output 
signal when the filter is applied on step and ramp like edges. 

Case 1: When is larger than the maximum grey value the region boundaries tend to smooth as non-linear Gaussian 

filter tends to work like a linear Gaussian filter. 

Case 2: When o-^ is large and is small the filter can sharpen ramp edges as shown in fig 5 below. 

It has been noticed that o-^ should be chosen greater than 1 to sharpen the ramp edge/. 

According to the authors a sequence of filtering operations are applied with varying parameters in order to 
achieve the desired smoothing effects. The first filtering step aims to reduce the contrast of the fine details in the 
image followed by the subsequent steps which does the same thing but tries to sharpen the edges of the coarser 
structures which is supposed to be blurred by the first filtering step. 


5 












I Procedia Computer Science 00 (2015) 7 -[ 4 ^ 


6 



original image Oj. = 0.5, a, = 128, rj = 1 


Figure 6. Smoothing the picture of a radio with five 5 filtering stages with o-x=0.5 and cr^=128. Figure reproduced from O. 


2.2.2. Effect of varying parameters 

In the first step cr^ should be very small so as to reduce the unwanted blurring of the coarser structures but at 
least larger than the size of the fine structures. The on the other hand determines the maximal contrast fine details 
approximately in order to be finally eliminated as a result. From the next filtering steps the cr^ needs to be increased 
and o-z to be decreased to sharpen the edges of the coarser structure. Experiments have proved that almost always best 
results are achieved if o^x is doubled and is halved before performing the next filtering step. 

2.2.3. Results 

Fig 6 below shows the result of filtering with a filter chain of five stages with initial parameters as (Tjc= 0.5 and 
(T^=128 respectively. 

This design introduces some fiexibility w.r.t its anisotropic counterpart as this Gaussian filtering can be consid¬ 
ered as an anisotropic diffusion process which can be stopped at certain times and can be restarted later with new 
parameters. 


f- E - p\\)ga,{fo{.q) - foip))-(f(q) - f(p))- (9) 

P qep 

where 

= 2 j - p\\)g<T,{fo{q) - foip))- ( 10 ) 

qep 

This operator smooths/ while retaining the edges present in/^. 

2.3. Nonlinear Gaussian filtering extension to bilateral filtering 
2.3.1. Gaussian smoothing to Bilateral filtering 

A Gaussian blur (a.k.a Gaussian smoothing) is the result of blurring an image by a Gaussian function typically 
used to reduce image noise and reduce detail. Each output image pixel value is a weighted sum of its neighbors in 
the input image. Mathematically, applying a Gaussian blur to an image is nothing but convolving the image with a 
Gaussian function. 

An image filtered by Gaussian Convolution is given by: 

G[I]p = Y,G.{\\p-q\\)Iq (11) 

qES 

where 

Gcr(x) = l/27icr^exp(-(x^ +y^)l2cr^) (12) 

denotes the 2D Gaussian kernel. It is a weighted average of the intensity of the adjacent positions with a weight 
decreasing with the spatial distance to the center position p. The weight for pixel q is defined by the Gaussian 


G^ = (\\p-q\\) 
6 


( 13 ) 
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a = 0.8 




a = 2 




a = 3 



Figure 7. Example of Gaussian linear filtering with varying cr. Top row shows the profile of a 2D Gaussian kernel and bottom row the result 
obtained by the corresponding 2D Gaussian convolution filtering. Edges are not preserved with high values of cr because averaging is performed 
over a much larger area. 


where cr is a parameter defining the neighborhood size. 


2.3.2. Drawbacks of Gaussian smoothing 

Suppose a bright pixel has a strong infiuence over an adjacent dark pixel although these two pixel intensity values 
are quite different. As a result, after convolving image edges are blurred because pixels across discontinuities are 
averaged together. Which concludes that action of the Gaussian convolution is independent of the image content but 
depends only their distance in the image. 

What is a bilateral filter? 

Before defining it at first the concept of domain and range filtering should be discussed, origin of which is linked 
to the nonlinear Gaussian filter as discussed above. Traditional filtering is domain filtering which focuses on the close¬ 
ness by weighing pixel values with coefficients which fall of with distance. Similarly range filtering averages image 
values with weight that decrease with intensity dissimilarity. It is this component which preserves the edges while 
averaging the homogeneous intensity regions. More likely to mention is that range filtering is non linear because their 
weights depend upon the neighborhood image intensity or color. Spatial locality is itself very promising and range 
filtering itself sometimes distorts an image’s color map. So when the domain as well as range filtering is combined 
together with both the advantages the duo is termed as bilateral filtering and is defined by, 

BF[I]p = -T - P\\)g(Trif(q) - f{p)) -if iq))-where (14) 

P qep 

= S<rs(\\q - p\\)g<rmq) - Kp))- (15) 

qep 

where Wp is the normalization factor, and denotes the standard deviation for the domain and range kernel, 
accordingly denotes the spatial Gaussian weighting whose value decreases with increase in pixel distance and 
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input spiatial kernel _f influence in the intensity weight 4" x g output 

■domain for the central pixel for the central pixel 


Figure 8. The bilateral filters. Each pixel is replaced by a weighted average of its neighbors. Each neighbor is weighted by a spatial component that 
penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures that only nearby similar 
pixels contribute to the final output. Courtesy of cni. 


Gcr^ which decreases with the influence of intensity difference i.e \f(q) - f(p)\ between pixels say p and q. Here we 
need to figure out the original framework from where the bilateral Alter equation has been derived. So we need to look 
back in equation and equation where see that the additional term i//(f(q) - f(p)) that has been added in equation 
8 to reduce the influence of f(q) as \f(q) - f(p)\ increases and later oni// = gg-^ is assigned. So this is actually the 
basis of the similarity function whose value is high in homogeneous regions within a boundary and drastically reduces 
across boundary on opposite sides of an edge and is analogous to the gg-^ in equation 13 which is nothing but the range 
kernel which is responsible for the sensitiveness to edges. 

2.3.3. Filter controlling parameters 

Unlike anisotropic diffusion Altering, the bilateral Alter is controlled by only two parameters to control its behavior 
namely the domain((T^) and range kernel((Tr) standard deviations respectively. 

Range kernelstd, dev.(crr). With the increase in cr^ the Alter gradually approximates towards the normal Gaussian 
convolution. As the range component increases it flattens and widens resulting in an increase in dynamic range of the 
image intensity so its very likely to penalize the identification of intensity gradient present at the edges. 

Domain kernel std. dev.(crs). The effect of variation of domain kernel std. dev. as already been shown in fig. 

It smooths the large features and blurs the detailed textures upon increasing. The response of the bilateral Alter is the 
multiplication of its constituent domain and range kernel respectively. For a large spatial Gaussian multiplied with 
narrow range Gaussian achieves limited smoothing in spite of the large spatial extent. The range weight enforces a 
strict preservation of the gradient contours. 

2.3.4. Iteration ejfect 

Bilateral Alter is inherently non-iterative in nature. However it can be iterated. Iteration leads to almost piecewise 
constant result as shown in fig. Likewise we try to relate the effect of iteration with the varying Alter parameters 
like o-y and cr^. It has been observed that since and are interlinked the effect of varying one doesn’t reflect much 
without the other parameter as shown in fig. 

2.3.5. Computational complexity 

There is no doubt that the bilateral Alter is a very cost effective algorithm to compute specially when the is large 
as it needs to compute a huge neighborhood of pixels, compute the weights of two parameters and cr^, followed by 
their product and an expensive normalization step. 

2.3.6. Decomposition into multiple components 

As shown in fig. [^the bilateral Alter can decompose an image into two parts namely the Altered image and the 
subtracted/residual component. The Altered image contains the large scale features, sacriflcing the minute textures 
respecting the strong edges. The remnant image left after subtracting the Altered image with the original one, contains 
only those portions the Alter removed, which can be treated further as either as flne textures or noise as the requirement 
demands. So we see that bilateral Alter provides an efhcient way to smooth an image while preserving the edges and 
has got various efhcient numerical schemes to realize the same together with its wide variety of application areas to 
be discussed in detail in the next subsequent sections. 
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Figure 9. The bilateral filters. Each pixel is replaced by a weighted average of its neighbors. Each neighbor is weighted by a spatial component that 
penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures that only nearby similar 
pixels contribute to the final output. Courtesy of Qo). 



1 iteration 2 iterations ^ iterations 


Figure 10. The bilateral filters. Each pixel is replaced by a weighted average of its neighbors. Each neighbor is weighted by a spatial component 
that penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures that only nearby similar 
pixels contribute to the final output. Courtesy of Go). 



1. Input .2. bilateral filtered 


3. Residual 


Figure 11. The bilateral filters. Each pixel is replaced by a weighted average of its neighbors. Each neighbor is weighted by a spatial component 
that penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures that only nearby similar 
pixels contribute to the final output. Courtesy of Qo). 
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(a) Extending all the regions (b) pixels within the neighborhood 

Figure 12. The bilateral filters. Each pixel in (a) is replaced by a weighted average of its neighbors as shown in (b). Each neighbor is weighted by 
a spatial component that penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures 
that only nearby similar pixels contribute to the final output. Courtesy of (H. 

3. Efficient ways to implement bilateral filtering 

Normal conventional bilateral filter as proposed by C.Tomasi and R.Manduchi Q is very computationally inten¬ 
sive especially when spatial kernel is large. Referring to equation 16 for each pixel p and q, whose loop is nested 
within p. So for a 10 megapixel photo there is 100 * 10^^ iterations involved which make the system very slow and 
would approximately take more than 10 minutes to process per image. 

To achieve the filtering there are several efficient approaches which reply on mathematical approximations yielding 
similar results with acceleration and accuracy that needs to be discussed in the following subsequent sections. 

- p\\)g<rxm - mum) (i6) 


3.1. Brute Force Approach 

The Brute Force approach follows the conventional way (direct implementation of the bilateral filter) as discussed 
with equations 16 to 18. Its computational complexity is (9(15 P) where \S \ the size of the spatial domain (cardinality of 
S) denoting the number of pixels. It’s clear that the quadratic complexity enhances the computational cost especially 
for large spatial domains. So there is a slight modification made with a view to reduce the computational complexity. 
The idea is to consider only the neighborhood pixels of the central pixel of interest p referring to equation 16 such 
that \\p - q\ \ < 2o-s considering the contribution of pixels outside the range of 2(T^ is negligible because of spatial 
kernel. The complexity has been decreased to (9(|51 crj) and is only effective for small kernels due to its dependency 
on o-g. But for this approach the computational complexity increases to a great extent. For a 8 megapixel image the 
number of iterations required is 64 * 10^^ which is obviously not a feasible solution (for the classical approach) for 
implementation in hardware as the computational complexity increases for processing nested loops of pixels. 

3.2. Jacobi algorithm approach 

In this section we need to discuss the ways to speed-up the bilateral filter and increase it smoothing effect together 
with its implementation procedure for piece-wise linear signals. Author here shown that the bilateral filter can be 
achieved from the Bayesian approach using a penalty function, and from it a single iteration of the Jacobi algorithm 
(also known as the diagonal normalized steepest descent DNSD) yields the bilateral filter. The corresponding penalty 
function is defined as 


N 

e{X} = \!2[X-Yf[X-Y]+XI2YjX-irxYW(Y,n)[X-irX\. (17) 

n=\ 

where X is the unknown signal and the result should be as close as possible to the measured signal Y. The matrix 
D stands for a one-sample shift to the right operation (161 . W is the normalized weight matrix. The bilateral filter can 
be accelerated in the following way: 
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A quadratic penalty function say: 


j 

e W = - P]X + Cj] (18) 

7=1 

the SD (steepest descent) reads 


" “Iz ^ - Pj] (19) 

7=1 

The author have shown a close resemblance of equation 19 with 20 by choosing , . from 20 as 2 1 = I, Qj^ = 

(/ - D^-^fW(Y,k - 1)(/ - D^-^) for k = 2,3,. ,N + 1 AND Pi = Y, Pi = P2 = P3 = .... = Pn+i = 0 for 

k = 2,3, . ,A+1. 

One effective way to accelerate the SD convergence is the Gauss-Siedel(GS) approach where the samples of Xi are 
sequentially computed from Xi[l] to Xi[L] (L scalar samples in vector Xi) and for the calculation of Xi[k], updated 
values of Xi are used instead of Xq values. This technique is known to be stable and converge to the same global 
minimum point of the penalty function given in 20. The GS intuition is to compute the output sample by sample and 
to utilize the already computed values whenever possible. 

While some describe this process in a more systematic way by decomposition of the Hessian to the upper-triangle, 
lower-triangle and diagonal parts. 

Alternatively the filter can be accelerated by slicing the gradient into several parts. From equation 21 it can be written 
for J iterations in the following form : 


Xi=Xo-ii[QiXo-Pi] 

Xi = Xi - ii[Q 2 Xi - P 2 ] 

X3=X2-li[Q3X2-P3] 


Xj = Xh - iilQjXli - Pj] (20) 

Proceeding in this way the final solution Xj tends to be closer to the global minimum point of the penalty function 
in 20. Finally a good de-noising has been achieved with reduced computational complexity but the computational 
analysis of the complexity has not been shown in the paper|[T6|. Well as of its implementation in hardware the itera¬ 
tive steps can be thought of to design as recursive procedure by loop unrolling in hardware. 


3.3. Layered Approach to approximate the filter 

Durand et. al |[T0l[T71 in order to accelerate the process of filtering, down-sampled the image 1 at lower resolution 
followed by computing the product with the range weight Gcr, and Gcr, the spatial kernel. The final layers are obtained 
by upsampling the convolution outputs followed by linearly interpolating between the two nearest layers as shown in 
step 3 of the pseudocode of the Algorithmic which dramatically accelerates the computation. 

As per Algorithm [C the downsampling process can be done in software domain itself. However as in step 1 the 
convolution operation of each layer with the spatial kernel can be accomplished in parallel with the specialized 
customized hardware filters meant for convolution operation as discussed in section 3.5 and 3.8 respectively which 
can accelerate the execution as the layers are independent to each other. Followed by the normalization operations 
where the pixels need to be divided, CORDIC divider can be applied thereafter which has got its inherent parallelism 
for the desired accelerated operation. 
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Algorithm 1 Layered Algorithm as proposed by Durand et. al tlOl 


Input: A 2D Image 1. A set of intensities {/q, h, •••in) is chosen and a set of layer |Lo, Li, .L„ | is computed as 

shown: = Gcr,(\ik - where / is a low resolution version of an image computed after downsampling. 

1. Each layer needs to be convolved with the spatial kernel followed by normalisation. 

Lk = {Gcr, 0 Lk) -i- {Go-, 0 Gcr,) where the denominator corresponds to the sum of the weights at each pixel and 
division is the per pixel division. 

2. The layers Lk is upsampled to get Lk. 

3. For each pixel p with intensity Ip, two closest intensity values 4i and ik 2 needs to be found and the output linear 
interpolation : 


BF[I,] ^ + 


ik2 Ip 




lk2-lk\ lk2-lk\ 

Return: Filtered image BF[Ip]. 


Bilateral 

Grid 



Gaussian blur 
grid values 

Slice: query grid 
with input image 



Filtered scanline 


(a) 



Downsampling in 2D Bilateral Grid 


(b) 


(c) 


Figure 13. The bilateral filters. Each pixel is replaced by a weighted average of its neighbors. Each neighbor is weighted by a spatial component 
that penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures that only nearby similar 
pixels contribute to the final output. Figure reproduced from J.Chen 


3.4. Bilateral Grid (3D kernel) approach 

J.Chen et. al ID evolved with a new idea of data structure to represent image data such that the kernel coefficient 
weights depends only on the distance between points. The new data structure is the bilateral grid (a volumetric data 
structure), enabling fast edge preserving image processing. Here a 2D image I is represented as a 3D grid T where 
the first two dimensions of the grid corresponds to the pixel position p and the third dimension the gray level pixel 
intensity Ip. The total process has been illustrated in Figj^and 14 with the corresponding algorithm in Algorithm 
The complexity of the algorithm is (9(151 + where |5| and \R\ are the size of the spatial and range domain 

respectively. There is a problem for color images as it becomes 5D grid and requires large amount of memory. 

Chen et, al ||8l have already implemented it in real time video rate using a dedicated graphics processor. Although it 
computationally very fast still we do not intent to use this procedure of filtering as it does not yield good denoised 
image quality in signal to noise ratio. 


3.5. Separable Filter Kernel Approach 

Pham et al. ID presented a separable implementation of the bilateral filter. The crux of the idea was that multi¬ 
dimensional kernel implementation is computationally expensive. So in order to approximate the 2D bilateral filter 
it is decomposed into two ID filters, filtering row followed by column or vice versa. Of course not all kernels are 
separable. Say at first the row filtering is applied and to the intermediate results the column filtering is applied to get 
the same 2D filtered output. The computational complexity decreases linearly from quadratic nature of brute force 
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Figure 14. The bilateral filters. Each pixel is replaced by a weighted average of its neighbors. Each neighbor is weighted by a spatial component 
that penalizes distant pixels and range component that penalizes pixels with a different intensity. Their combination ensures that only nearby similar 
pixels contribute to the final output. Figure reproduced from J.Chen (3. 
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Algorithm 2 Bilateral Grid algorithm 

Input: A 2D Image I. A grid Y has been build such that Y : SXR containing homogeneous values shown in 

Figli^a. 

r{px, Py, r) = (Kpx, Py), 1 ) if r = I{px, Py) 

=(0,0) otherwise. 

1. Each layer needs to be convolved with the spatial kernel followed by normalization. 

Lk = (Gcr, 0 Lk) ^ (Gcr, 0 Gcr,) where the denominator corresponds to the sum of the weights at each pixel and 
division is the per pixel division shown in figf^c,d. 

2. Downsample Y to get f. in fig[^b and|l4(b respectively. 

3. Perform a Gaussian convolution off for each component as shown: 

GC[Y](p^, Py, r) = Gcr^cTr ® f (Px, Py. r), 

where Gcr,,crr is a 3D Gaussian, and are parameters along the spatial and range dimensions respectively. 

4. Upsample GC[f] to get f. 

Return Result: For a pixel p with initial intensity Ip, (wl, w) the value at position (p^, Py, Ip) in f is denoted. The 
result of the bilateral filter is BF[I]p ^ wl/w. 


approach and leads to 0(\S | cr^) as it neighborhood computation is single dimensional, where \S | the size of the spatial 
domain (cardinality of S) denoting the number of pixels and cTs defines the neighborhood size of the spatial domain. 
Separable kernel computes the axis-aligned separable approximation of the bilateral filter which is suitable for ho¬ 
mogeneous areas and symmetrical edges but performs poorly in regions of complex textures. There are methods to 
overcome such problems but it comes with an additional computation overhead. We have successfully implemented 
the separable kernel based concept of 0 applied on the trigonometric based kernel mil on FPGA and obtained a 
good acceleration as shown in llT9l . 

3.6. Local Histogram Approach 

Bein Weiss CD computed the bilateral filter from the histogram of image. Considering the spatial weight as a 
square box function the bilateral filter can be represented as containing only the range kernel component with the 
domain kernel as a constant entity. Their idea is as follows: Considering two adjacent pixels p\ and p2 and their 
corresponding neighborhoods Ng-XP^) ^nd Ncr^P'^)- Based on these scenario the histogram of the neighborhood 
No-Xp^) is computed keeping an eye with the similarity of the histogram of the other neighborhood. Once it is known 
the bilateral filter can be computed as shown: 


BF[I]pi = 1/Wpi Yj 

Tf\lpX - Ip2\)Ipl 

(21) 

p2eN^,(pl) 



II 

GaX\lpl - Ip2\) 

(22) 


p2€N^,(pl) 


This implementation if followed by iteration the filter thrice to remove the strange band artifacts near edges of 
strong intensity gradients which remains after filtering. The complexity of the algorithm is of the order of 0(151 logo's). 
It can handle filter kernels of any size in short time but has a setback w.r.t processing color images independently for 
each separate channel. It can exploit vector instructions of CPU and as a result parallel operations of the independent 
instruction can be implemented concurrently with full acceleration. Interested readers need to refer to the paper na 
for the detail explanation of the algorithm. 

3.7. Constant time bilateral filter implementation using polynomial range kernels. 

Eater on Porikli Q in 2008 showed that Gaussian range and arbitrary spatial bilateral filters can be realized by 
Taylor series as linear filter decomposition with a noticeably good filter response. Considering arbitrary spatial filters 
a polynomial range has been defined as: 
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g(p+k) = [i-(i(p)-i(p+k)fr ( 23 ) 

where n is the order of the polynomial. Considering bilateral filter of the form 

yb(p) = h~^Yj -^ip + k)) (24) 

keS 

where g{p, k) and f{k) denotes the range and domain kernels respectively and the normalising term k^ = Yj f(k)g(I(p)- 
I(p + k)). Now for ^ = 1, in Eqn 25 and substituting the value in Egn 26 the polynomial range arbitrary spatial filter 
is obtained as 

yb(p) = kC [E mnp + k)- l\p) Yj f(k)I(P + k) + 2 I(P) Y mi\p + k)-Y nk)l\p + k)] (25) 

The corresponding linear filter responses are yi = 2 f(k)I(p + k), y2 = 1] f(k)I^(p + k), etc. Thereby substituting 
the above values in Eqn. 27 it can be rewritten as 


yb = kb ' [(1 - I^)yi + 2Iy2 - ys] (26) 

where = I(p), = I(p)I(p), etc. Accordingly the normalization term is kt, = I - 2Iyi - y 2 and the spatial 

filter / is not constrained here resulting for any desired filter to get chosen. Eor {n = 2) Eqn. 25 can be written from 
Eqn. 26 as 


yb = kb-^ [(1 - 2f + P)yi + 4(1 - P)y 2 + (6f - 2 )j 3 - 4Iy^ + y^] (27) 

where kb~^ = l-2f + n + 4(7 - fi)yi + (6f - 2)>>2 - 41y2 + >’4 

So its clear from equations 28 and 29 that the bilateral filter has been achieved with polynomial range terms 
in terms of the spatial filter without approximation. The author had presented one more way namely the Gaussian 
range. Additional smoothness to realize the bilateral filter by exploiting the property of Gaussian filters to be easily 
differentiable and can be expressed in terms of linear transforms. The range component is given by 

exp(-a[I{k + p)-I(p)f). (28) 

The equation can be rewritten as 

exp(-afi(p))exp(-a[fi(p k) - 2I(p)I(p + k)]). (29) 

Thereby applying Taylor expansion to Eqn. 31 the bilateral filter expansion upto second order and third order 
derivatives can be obtained as shown: 

yb ~ kb~^lyi + 2(27^2 + (x{2afi - l)y 3 - 2a^Iy/[ + 0.5Qf^y5] (30) 

yb ~ kb~^[y\ + 2aly2 + (2a^f - a)y 2 - 2a^(l - ^al^)y4 + ^^(O.S - 2fa)y5 + a^Iye - ia^l6)yj] (31) 

Additively the normalizing terms are similar to some extent containing same terms. So the bilateral filter has been 
realized using the weighted sum of the spatial filtered responses of the powers of the original image. Their complexity 
does not scale with the shape or size of the filter and are referred to as 0(1) algorithm as its complexity. 

It is not a good idea to implement this algorithm in hardware as it would increase the complexity for developing the 
power series as exponential operation execution is not feasible in hardware. 
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3.8. Constant time bilateral filter implementation using trigonometric range kernels. 

Here the authors Kunal et. al CD have used the class of trigonometric kernels which turned out to be sufficiently 
rich, allowing for the approximation of the standard Gaussian bilateral filter. This has been done by generalizing 
the idea of the using polynomial kernels (71. It has been observed that for a fixed number of terms the quality of 
approximation achieved using trigonometric kernels is much superior than obtained using polynomial kernels. This 
was done by approximating the Gaussian range kernel of the bilateral filter using raised cosines. They observed 
that trigonometric functions share a common property of polynomials which allows one to linearize the otherwise 
non-linear bilateral filter. The trigonometric functions share a common property of polynomial which allows one to 
linearise the otherwise nonlinear bilateral filter. Recalling the original equation of the bilateral filter in Eqn 16, it is 
seen that the nonlinear property of the bilateral filter mainly arises due to the range kernel component. If the range 
part is kept constant and is replaced by a unity value then the filter can be implemented in constant-time, irrespective 
of the shape or size of the filter. And it is only possible if the range component can be written in the form of a cosine 
term say 


^cr,(s) = cosiys) (32) 

where {-T < s < T) and T is the dynamic range of a grayscale image and is equal to 255 and y = njlT. The 
authors have also shown that this idea can also be applied to more general trigonometric ideas. Interested readers can 
view the details. Now in order to penalize large difference in intensity than small ones it is needed that the family of 
cosines should be symmetric, non-negative and monotonic. These properties are all enjoyed by the family of raised 
cosines of the form 


gcr^(5’) = [c6>5'(75’)]^ where (-T < s < T) (33) 

Now as the value of the raised power N in Eqn. 35 gets large the go-X^) quickly gains Gaussian like behavioral 
curve over half the period [-tt, tt] as N increases. This property speeds up the rate of convergence than that of the 
Taylor polynomials Cl used to approximate the Gaussian range kernel by suitably scaling the Gaussian kernel. This 
approach greatly reduces the computational complexity by 0(1) and can also be implemented in parallel further ac¬ 
celerating its speed producing artifact-free results. We like to refer to his article ifTlI for the detail of the algorithm. 
We have further improved it w.r.t speed of execution by realizing the filter on an EPGA CS). We have done an EPGA 
based implementation of an edge preserving fast bilateral filter on a hardware software co-design environment of this 
most recent algorithm preserving the boundaries, spikes and canyons in presence of noise. Trigonometric terms (like 
sine and cosine as in equation 34 and 35) implementation meant for trigonometric range kernel execution for realizing 
the fast bilateral filter has been done using the customized circuit (with CORDIC processor) for robust angle computa¬ 
tion. Customized CORDIC processor has been designed to compute sine and cosine data with a provision to compute 
large angles more than the inherent tt and 2 * tt to accommodate/map within it. Eollowed by the four stage parallel 
pipelined architecture greatly improves the speed of operation. Moreover, our separable kernel implementation of 
the filtering hardware increases the speed of execution by almost five times than the traditional convolution filtering, 
while utilizing less hardware resource csi. Author in CD have shown that by exploiting the central property of 
trigonometric function aka shiftability which assists in reducing the redundant computations inherent in the filtering 
operations. It has been shown that using the shiftable kernels certain complex filtering can be reduced to simply that 
of computing the moving sum of a stack of images where each image in the stack is obtained through an elementary 
pointwise transform of the input image adding two significant advantages. Eirstly fast recursive algorithms can be used 
for computing the moving sum and, secondly, parallel computation can be used to further speed up the computation 

CUIISl. 

Some of the benchmark implementation outcomes has been presented in Eigures[^and[T^respectively. Erom Eig. 
p3]we get that the output of Tomasi's method produces a particular PS NR value of the denoised output but consumes 
a huge time. The outcome of Chen’s method in (d) takes the least time of all but with a cost of a decreased PS NR 
with a slight degraded image quality, (e) shows better results w.r.t both the PS NR and time of execution. Also the 
quality of approximation achieved using trigonometric kernels is much better than obtained using polynomial for a 
fixed number of terms. Depending on the availability of the resources, the Bilateral grid approach of Chen (Sj executes 
well for high resolution real time videos on graphics hardware. The local histogram approach of Weiss ifT^ is better 
suited for filter kernels of arbitrary size. Einally as discussed above a step ahead to reduce the complexity to 0(1) 
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(a) (b) (c) (d) (e) 


Figure 15. (a) denotes the original image, (b) denotes the noisy image with added noise of cr = 0.15, PS NR = 20.42db, (c) shows the denoised 
image after applying bilateral filtering with Tomasi's method ( 3 , PS NR = 21.01db, time taken to execute is 1.28 second, (d) shows the denoised 
image after applying bilateral filtering with Chen’s (3, PS NR = 20.40db, time taken to execute is 0.013 second, (e) shows the denoised image 
after applying bilateral filtering with Kunal's method (m , PS NR = 2l.l4db, time taken to execute is 0.75 second. 



Figure 16. The output of two different bilateral filters. The left hand side figQ^a shows the output of method Q and right hand side fig[T^b is of 
m Note that the strange artifacts in the L.H.S figure particularly around the right eye (see zoomed insets) is due to the approximation caused due 
to polynomial unlike the one in R.H.S. Figure reproduced from K.N Chaudhury et. al E]. 


Taylor series polynomial followed by a more advanced version the trigonometric range kernel has been deduced so 
far. Tablesummarizes the complexity order of each of the previously discussed techniques. 


Table 1. Complexity analysis report 

Algorithm Complexity 

Brut Force Approach (section 3.1) 0(\S\^) 

Layered Approach (section 3.3) 0(151 + ^ 

Bilateral grid(3D kernel) approach (section 3.4) 0(|51 + ^ ^) 

Separable filter kernel approach (section 3.5) 0(|51 cr^) 

Local Histogram Approach (section 3.6) 0(|51 logo's) 

Constant time polynomial range approach (section 3.7) 0(1) 

Constant time trigonometric range approach (section 3.8) 0(1) 


4. Applications of edge preserving filters 

This section discusses the the uses of the edge preserving filters (anisotropic and bilateral) for a wide range of 
applications. 

4.1 . Contrast Enhancement. 

There lies some problems to eliminate the noise from the edges in case of nonlinear isotropic diffusion filter. 
Anisotropic diffusion do solve the problem as it considers both the modulus of the edge detector |v£’o-| as well as its 
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Figure 17. Anisotropic diffusion filtering applied on a Gaussian type function with parameters cr = 2, A = 3.6 and Q = (0,256)^. Top to bottom all 
the figures are denoted with time beginning from t=0, 125, 615, 3125, 15625, 78125. Figure reproduced from B.Stuttgart et. al do). 


direction. The authors (201 constructed an orthonormal system of eigenvectors vl, v2 of the diffusion tensor D such 
that 


VilIV^^, V2-LV^^ 


(34) 


It is needed to be smoothed along the edge than across it. As a result Weickert BSll choosed two eigenvalues Ai 
and A 2 such that 


Ai(vE^) := g(\vE^\^) (35) 

A2(vE^) := 1 (36) 

The goal of the filter decides the choice of the eigenvalues Ai and A 2 . Smoothing within each region and inhibition 
of diffusion across edges can be obtained by reducing the diffusivity Ai perpendicular to edges. Increasing it enhances 
the contrast. It depicts the changes in the Gaussian like function with temporal evolution. This edge location remains 
stable over some time till the process of contrast enhancement is concluded the steepness of edge decreases and the 
image converges towards a constant image. 

4.2. Coherence enhancement 

It begins with the wish to process one-dimensional features like line-like features |[2Ql . Cottet et. al 1^ with their 
diffusion tensor D, eigenvectors as in Equation 36, with the eigenvalues as shown: 


l+(|v£^|/cr)2 

18 


A2(VE^) = 


Ti(v£^) = 0 
(r] > 0). 


(37) 

(38) 
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{a) (b) 


Figure 18. Anisotropic diffusion filtering applied on a fingerprint image, (a) image shows the distorted noisy input image while (b) the coherence 
enhanced image applied with parameters, cr = 0.5,p = 4, t = 20. Figure reproduced from B.Stuttgart et. al (20). 


which denotes the diffusion process diffusing perpendicular to Vw^-. As cr ^ 0, Vu becomes an eigenvector of 
D{jp{yua-)) with the corresponding eigenvalue diminishing to 0, halting the process completely. Now by substituting 
v£’o- by a more robust local descriptor one can improve the parallel line like textures thereby enhancing the coherence¬ 
enhancing anisotropic diffusion. Figureshows the results of the AD filtering applied on a fingerprint image. 

4.3. Denoising 

Both anisotropic as well as bilateral filtering are used for denoising purpose. Here in this section we try to focus 
more on the use of bilateral filter to remove noise from digital images. 

There are several sources of noise corrupting a digital image like dark current noise ll^ due to thermally generated 
electrons at the sensor site, shot noise arising out of quantum uncertainty in photo-electron generation. Amplifier noise 
and quantization noise occur during the conversion of the number of electrons generated to pixel intensities (221. The 
basic purpose of the bilateral filter is image denoising and has a wide range of applications like image restoration, 
medical imaging, object tracking etc. to name a few. Bilateral filtering is very simple, easy to understand and its 
recent implementations are very fast. It preserves object contours (edge preserving) and generates sharp results. It 
smooths the homogeneous intensity regions while respecting the edges and in this way it reduces the noise injected in 
images but sometimes the edge preservation is not accurate and introduces undesirable staircase ejfect and cartoon¬ 
like features. It is undoubtedly a good solution but not for all kinds of application. The application areas where it has 
been applied are : 

4.3.1. Normal denoising 

As we know that bilateral filtering is used for denoising purposes. We like to validate it via some experimental 
results as shown in Figure [T^ and with the parameters correspondingly as shown. Figure [T^ shows the original 
image with its corresponding noisy one after adding white Gaussian noise. The right hand side contains its filtered 
version with zoomed insets surrounding it to provide a clear vision. Figurej^shows the filtered version applying var¬ 
ious benchmark methods starting from normal classical implementation as in a, h shows the 3D grid implementation 
and c shows the trigonometric kernel implementation ifTTIl . The corresponding denoised images are also visible and 
they are used as per the demand in applications. 

The experiment has been implemented 10 times in software and the corresponding mean squared error (MSE) 
obtained has been averaged by 10 to get the averaged MS Eav, which is used to calculate the PS NR of the software. 

4.3.2. Video implementation 

T.Pham m, j. Chen (81 Bennet et. al (23ll shows that bilateral filter can be implemented for real time videos. T. 
Pham have shown that better image quality and and compression efficiency can be achieved if the original video is 
preprocessed with the separable filter kernel computing at a fraction of executing time unlike the traditional one. Chen 
et. al have shown the use of a new data structure the bilateral grid approach already discussed in section 3.4. Using it 
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Figure 19. Filter output for image size of 150 * 150 for the additive Gaussian noise. Filter settings crs=20, crr=50 and cr=12 for the additive 
Gaussian noise. Figure reproduced from (m 




Figure 20. PSNR values of the filtered images for the different Gaussian Bilateral filters on the grayscale image of Einstein of size 150 * 150 with 
filter settings 0 -^= 16, cr^= 0.1 and cr= 0.15 for the additive Gaussian noise. Figure reproduced from m. 
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(a) without prepmcessing 


(b) with bilateral pre-filtering 


Figure 21. Frame 31 of a Foreman sequence compressed at 150 Kbits/s where the video without preprocessing on the left in (a) is more blocky 
than its counterpart shown in (b). The process is described in section 3.5. Figure reproduced from (3. 



(a) (b) (c) (d) 


Figure 22. (a) Input frame, (b) Histogram stretched version, (c) red and green shows the number of temporal and spatial pixels integration, (d) 
Output of Bennett et. al Ea. It describes the process to combine spatial and temporal bilateral filtering to achieve high-quality video denoising 
and exposure correction. Figure reproduced from (23). 


edge aware image manipulations such as local tone mapping on high resolution images can be done at real time frame 
rate with parallelising the algorithmic instructions on modern GPUs. 

4.3.3. Medical Image Processing 

In medical image processing Wong et. al El and D. Bhonsle et. al ES removes the noises corrupting the 
images mainly in presence of white Gaussian noise. Although it doesn’t work good with salt and pepper noise. Wong 
improved the structure preservation abilities of the bilateral filter by an additional weight component depending on 
the local shape and orientation known as trilateral filtering. 

4.3.4. Orientation smoothing 

D.Chen et. al (261 introduced first the concept of orientation bilateral filter which combines the concept of 
average(arithmetic mean value) orientation and the bilateral filter to maintain micro structural features and to smooth 
orientation noise which needs a fast computation of the disorientation using quaternion equation by Hamiltonian as 
follows: 


Q = cos(6/2), risin(6/2), r2sin(6l2), r3sin(6l2) (39) 

where r = (ri, r 2 , r 3 ) denotes the is the axis and 0 is the angle of rotation of the axis which completes the orientation 
specification. The average orientation based on the orientation filter is given by the weighted function: 
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Q 


f 



( 2 ) 


Figure 23. (1) Hexagonal sampling grids meant for the bilateral filtering process of a hexagonal image structure, (2) (a) orientation maps for 
the cross-section of electro-deposited copper using inverse pole figure color coding with the reference direction in the normal direction (ND): (b) 
without a filtering process, (c) using a hexagonal sampling grid of 7 sampling points, (d) using a hexagonal sampling grid of 19 sampling points, 
(e) using a hexagonal sampling grid of 37 sampling points, and (f) the inverse pole figure color coding. Figure reproduced from ED. 


Q(xo,yo) 


Zto Q(Xi,yi) * W(Xi,yi) 
Z":o W(Xi,yi) 


(40) 


where the weight W(xi,yi) is the product of two components namely the domain and range part respectively 
as given by Ws(xi,yi) and Wpixi^yi) which are obtained by methods described in section 2.3.1 of this paper. The 
intensity difference of the range component here denotes the misorientation between the sampling point and 

the center point (vo,yo)- the smoothing process is continued as per the equation 42 and the output is as shown in 
Figure]^ Through this the authors have tried to improve the angular precision of orientation maps for deposited and 
deformed structures of pure copper obtained from electron backscattered diffraction (EBSD) measurements used in 
the quantitative characterization of crystallographic micro structures. 

Paris et. al used also applied the same filter to smooth the 2D orientation field from optical measurements for hairstyle 
modeling WT\ as shown in Figure [24| 


4.4. Contrast management 

Bilateral filtering can be applied for contrast management applications like grouping large and small scale com¬ 
ponents of an image (detail enhancement) or reduction, texture and illuminating separation, tone mapping and man¬ 
agement and generation of high dynamic range imaging na. 


4.4.1. Decomposition for detail enhancement and Image fusion 

For detail enhancement Fattal et. al (23 used a new image-based technique for enhancing the shape and surface 
details of an object. They computed a multiscale decomposition based on the bilateral filter for a set of photographs 
taken from a fixed viewpoint, but under varying lighting conditions as an input. For a basic multiscale bilateral de¬ 
composition its purpose is to build a series of filtered images V that preserve the strongest edges in I while smoothing 
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(a) Ztiom uti bpul image 


(b) Oricneatiuti^ before bilateral fiktring (c) Orienialioiiia after bilalera! filleriiig 


Figure 24. Hair orientation measurements with applied bilateral filter to a complex plane. Figure reproduced from Paris et. al im. 


small changes in intensity. At the finest scale 7 = 0 ,/® = / is set and then iteratively the bilateral filter is applied to 
compute: 



(41) 

qeQ. 


k = Y8crs.jm\)8<rr.jilU,-4) 

(42) 


qeO. 

where p is the pixel coordinate, cr^j and cr^j are the widths of the spatial and range Gaussians respectively and q 
is an offset relative to p that runs across the support of the spatial Gaussian. Finally an enhanced image combining 
detail information at each scale across all the input images needs to be reconstructed. It is simple to implement, 
fast 0(N^logN) and accurate. Results shown in Figure ^ It is to be noted that the bilateral filter is a good choice 
for multiscale implementation as it avoids the halo artifacts commonly associated with the classical Laplacian image 
pyramid. 

Image fusion is again categorized into two applications namely (a) flash photography and (b) fusion of multispec- 
tral images. 

Flash photography: Visually compelling images can be generated in low lighting conditions by combining flash and 
no-flash photographs. Elmar et. al 1281 used the bilateral filter to decompose the image into detail and large scale. 
The finer textures(detailed feature) extracted from the flash photography is combined with the large scale components 
of no-flash image. 

Multispectral fusion: Here both the data from infrared spectrum and standard RGB data is used to denoise video 
frames in low light conditions as per Bennett et al. ll2^ . A modified bilateral filter (dual bilateral filter) is applied in 
such condition where the range weight is constituted by a combination of both the visible spectrum and the infrared 
spectrum given by: 


BFARGB]p = HWp Y, GaXWP - - /?G5,||).G<,„(||//?p - IR,\\)RGB, (43) 

qeS 

where RGBp is a 3-vector representing RGB component at pixel p and IRp the intensity of the infrared spectrum 
at the same pixel p. 




rnf 


peQ. 


(44) 


where / and nf denotes the flash and no-flash pictures, and f(g - s) denotes the domain kernel and g(Ip - l{) 
denotes the range kernel. Here the edge preserving term although corrupted in the no-flash image is preserved. 
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no-flash flash result 

Figure 25. Figure showing the no-flash,flash and the resultant output of the previous duo. Figure reproduced from Elmar et. al m. 



Figure 26. Usage of bilateral filter to create multiscale decomposition of images revealing minute details after decomposition and combining 
shading information across all the images. Figure reproduced from Fattal et. al Ell 


4.4.2. Large and small scale component separation 

It has been observed that for large range kernel standard deviation cr^, the bilateral filter eliminates the fine texture 
variations of refiectance while maintaining large discontinuities of illumination changes 1^ . The crux of the idea is 
illumination changes occur at a large scale than texture patterns. 

4.4.3. Tone Mapping and management 

Durand et. al cni in the year 2001 demonstrated the significant property of bilateral filter to seclude small-scale 
signal variations, can be utilized to design a tone mapping process which can map intensities of high dynamic range 
image to a low dynamic range display. They overcame the limitations of the existing classical approach like gamma 
reduction, uniform scaling etc. which wipes out the finer details and texture depths of displays. Here the bilateral filter 
is applied on the log intensities of the high dynamic range images, followed by evenly scaling down the result and 
summing up filtered residual, resulting in a more visually pleasing display (TOl . The texture of the image is defined 
by: 


^ y G.Mp - q\\)G^X\logIp - %4|) \H\, (45) 

qeS 

where the texture can be varied by adjusting the small scale component visible in the image. For contrast man¬ 
agement purpose the bilateral filter is usually applied to the log of the source image as the response of human visual 
sense is multiplicative. Also logarithm operation helps the standard deviation of range component to act uniformly 
across different levels of intensity. However using a cubic root which is based on luminance channel of the CIE-Lab 
color space, handles more efficiently the multiplicative process due to its computational simplicity. 
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Figure 27. Usage of bilateral filter to create a tone mapped image. The six individual images in (a) is used to generate the image as shown in the 
right hand side (b). It is to be noted that the fine details of interior of the church is not visible due to underexposed and overexposed regions in light, 
and grabbing all the information from the six images from (a), figure (b) is generated where almost all of the fine details (depth) are visible. Figure 
reproduced from wiki. (30). 


4.5. Depth Reconstruction 

Bilateral filter improves in depth reconstruction, i.e recovery of pixel depth value from corresponding pixels point 
pairs of different views. Q.Yang et. al ED constructed a metric called cost volume, which is a 3D volume of depth 
probability, provided by the probabilistic distribution of depth of input range map image. Then they applied iteratively 
the bilateral filtering on the cost volume to generate visually compelling high resolution range images with accurate 
depth estimate exhibiting a very effective but simple formulation. Yang 1^ again in the year 2009 showed the usage 
of bilateral filter in stereo reconstruction. Among several strategies for depth reconstruction bilateral aggregation 
turned out to be the most successful one. 

4.6. Feature preserving mesh smoothing 

Thouis et al. ll33l showed a noniterative and robust mesh smoothing technique to remove noise while maintaining 
the feature using a robust statistical approach and local first order predictor on surface. Predictors have been used 
based on triangles of mesh on natural tangent planes to the surface. The readers should a the Figure. for clear 
understanding. A new predicted point p is being generated based on the predictor from its nearby spatial 

triangle. The spatial weight / depends on the distance ||/7 - c^|| between p and the centroid Cq of q. The range 
component g depends on the ||n^(p) ~ p\\ between the prediction and the original source point p. Therefore the 
prediction p on a surface is given by: 

p - ;’||)g(||n^(;’) - p||), (46) 

qes 

where the normalizing factor k is : 


k{p) = aJi^Cq - J3||)g(||n^(p) - (47) 

and Oq is the area of the triangles used to weight for the variations of sampling rate of surface. 

In Figure the normals as shown are more prone to noise than vertices as they are first-order properties of the mesh 
as shown in Figure [2^. However the estimation can be increased with mollification |[3^[35]i . which is achieved by 
smoothing the normals without changing the positions of the vertices during mollification to preserve the features as 
few comers might get improperly smoothed by mollification near corners. The results are shown in Figure 

There are other applications of bilateral filtering which includes video stylization (361, demosaicking (381, optical 
flow Ell etc. 


25 



I Procedia Computer Science 00 (2015) 7-p^ 


26 


0 riqfp] 



Figure 28. (a) A point p is predicted as Wq{p) based on the surface at q is the projection of p to the plane tangent to the surface at q. As we 
see the predictions are far away constructed from points across sharp features, (b) Noisy normals results in poor prediction (c) modification after 
mollification. Figure reproduced from Thouis et al. oa 



(a) (b) (c) (d) (a) 


Figure 29. (a) Original image with mesh, (b) Isotropic smoothing on. (c) approach without mollification, (d) approach with no influence weight 
g. (e) complete approach of Thouis oa. Figure reproduced from Thouis et al. Ea 


5. Connection among various edge preserving filters 

There exists a number of existing edge preserving smoothing filters having various contributions in graphics and 
image processing and among them we find bilateral filter and anisotropic diffusion as the most dominating one and 
will show their relationships with one another and with other methods as they give similar results. 

The bilateral filter acts as a bridge between anisotropic diffusion filtering and other filtering methods which is 
shown pictorially as in Figure]^ 

5.1. Yaroslavsky filter with bilateral filter 
Let us refer to an ingeneral equation as 


v(/) = u{i) + ^(0 (48) 

where v(/) is the noisy output image, u(i) is the pure ideal image and n(i) is the additive noise. Yaroslavsky filter 
averages the pixels with a similar grey level value and belonging to the spatial neighborhood Bp(x) given by 

1 r \u(y)-u{x)\^ 

YFh,pu(x) = —— I u(y)e dy, (49) 

Lw JBpix) 

p _ \u(y)-u(x)\^ 

where x6kl,C(x) = L . .e Jy is a normalizing factor and /z is a filtering parameter. Later on in 1995 

^ tip } 

and in 1998 SUSAN filter (S) and the bilateral filter ID got introduced which instead of considering a fixed spatial 
neighborhood Bp{x), weigh the distance to the reference pixel v. given by 

1 r \u{y)-u{x)\^ 

BFh,pU(x) = u(y)e d e dy, 

C(x) Jn 

26 


( 50 ) 















I Procedia Computer Science 00 (2015) 7-p^ 


27 



Figure 30. Relation among various edge preserving filters. 


As per Baudes in 2005 11521 there is no as such difference between bilateral filter(BF) and Yaroslavsky filter(YF). 
Bilateral filter with a square box window is simply an approximation of the Yaroslavsky Filter which restricts the 
integral to a spatial neighborhood independent of the position y. In other words Yaroslavsky filter is a special case of 
Bilateral filter with a step function as the spatial weight. As a result we don’t find the term |y - xf' in the Yaroslavsky 
filter in equation [^unlike in[^ If the gray level intensity difference is larger than h, then both the filters compute 
the averages of the pixels belonging to the same region of the reference pixel thereby behaves as an edge preserving 
smoothing filter. 


5.2. Relation of bilateral filter with anisotropic dijfusion filter 

There exists a wise link between anisotropic diffusion filtering and bilateral filter and here we would like to focus 
on the discrete domain computation. Fortunately it has been observed by D.Barash (401 that adaptive smoothing 
serves as a link between anisotropic diffusion and bilateral filtering. 

Anisotropic diffusion and adaptive filtering: 

Give an image Pfix), where ~x = (xl,x2) are the spatial coordinates, and the iteration of the adaptive filtering 
yields: 






+/,x2 + j)w‘ 


where the convolution mask is defined as : 


(51) 


\dfxl, v2)f 

w\xl, x2) = expi-- - ——-) (52) 

k the variance of Gaussian kernel, where d\xl, x2) is given by the magnitude of the gradient computed as: 

d\x\,x2)= (53) 

where 


{Gxi^Gxi) = ( 


dl\x\,x2) dl\x\,x2) 
dxl ’ dx2 


(54) 
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which bears similarity between the convolution kernel of the diffusion coefficient of the anisotropic diffusion. 
Considering one dimensional signal of I^x) and showing the averaging process as follows: 

/'+I(x) = Cil‘ix - 1) + C 2 l\x) + C 2 l\x + 1) (55) 

where 

Cl + C 2 + C 3 = 1 (56) 

which can also be written as, 

I‘^\x) - I\x) = Ci(/'(x - 1) - I\x)) + C 3 (/'(x + 1) - l\x)) (57) 

replacing ci = C 3 the above equation reduces to: 

P^\x) - I\x) = Ci(/'(x - 1) - 2l\x) + I‘(x + 1)) (58) 

which can be written as a discrete approximation of the linear diffusion equation as: 

(59) 

ot 

Considering a second case when the weights are space-dependent: 

r^\x) - I\x) = c\x + l)(/'(x + 1 ) - I\x)) - c\x - \){l\x) - I\x - 1 )) (60) 

which finds similarity with the implementation with the anisotropic diffusion equation by Perona and Malick (2) 

^ = v(c(xl, x2) V /) (61) 

Ot 

where 

c(Xi,X 2 ) = g(||V/(Xi,X 2 )||) (62) 

and ||V/|| is the magnitude of the gradient and g(||V/||) is an edge stopping function satisfying every parameters. 
Thus a link between anisotropic diffusion(60) and adaptive filtering(50) has been established. 


Bilateral filter and Adaptive filtering. 

The discrete version of the bilateral filter is given below: 


Ills IlLs^‘ 


(63) 


where the weight is given by: 


w (V, = exp (— —-r )exp( --) 


2crl 


2<tI 


(64) 


where S is the kernel size of the filter and the generalized intensity consisting of domain and range components is 
given by: 


/ = 



( 65 ) 
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Figure 31. (a) A selected window and its neighborhood meant for histogram calculation, (b) As per the figure each pixel tends to shift towards the 
maximum of the local mode, (c) Curve shows the smoothing effect of the local histogram upon change in the range parameter, concept of local 
mode filtering reproduced from EH- 


exp(- 


w‘Cx) = exp{-- 


1(0 ^ 


Cr cr d 


HO-H^) 

( /(x) ~x 1 


O'R O' J) 


) 


exp{- 




O-R 


O-D 


^ H(HO- 

= exp(--( -+-)) 

2 crR cr D 

= exp( -— )exp(, -—^-) 


2cri 


2crl 


( 66 ) 


Thus it can be seen from equation 65 that it contains two components namely the domain and range kernel weights 
respectively which proves the link between bilateral filtering and adaptive filtering. Thus it can be concluded that since 
there is a link between bilateral filter and anisotropic diffusion filter via adaptive smoothing filter so knowing one the 
other can be derived. 


5.3. Bilateral filter with local mode filtering and robust mode filtering 

Sylvain et al. nil have observed that starting from local mode filtering, the bilateral filtering can be derived. It can 
also be achieved while trying to find a solution of a robust minimization problem. Here we like to summarize it. The 
local mode is basically chosen unlike the global mode(shown in Figure |3T]:.) to preserve the fine details of an image. 

Given a grayscale image I the local mode filtering can be defined as: 

= (67) 

qeS 

such that / : Q ^ d is the Dirac delta function. Since the local mode filtering also depends upon the smoothness 
parameter of the histogram, so it can be defined as. 


H\i, <Tr) = 0 GHl) = Yj ^<^r(h - 0, 

qeS 


( 68 ) 
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where cr^ is the intensity range within which the smoothing is done. Now to define a histogram locally i.e around 
a given position p a spatial Gaussian kernel needs to be applied around the point p to get: 

H\p, i, CT„ cr,) = Yj G^siWP - q\\)G^P, - 0, (69) 

qeS 

where as denotes the spatial neighborhood centred around pixel position p. 

As all the pixels within the local neighborhood tends to move towards the maximum local range centered at Ip so 
it is required to equate the equation to zero as shown: 

dH^ 

—(p, i, (Tr, O-s) \i=I^= 0. (70) 

The equation 69 equates to be: 


J](4 - i)G^x\\p - q\\)GcrXi, - 0 = 0, (71) 

qeS 

where Ip needs to verify the following implicit equation, equating Ip = i where i is: 

^ i:,esG<rMP-q\\)G<rP,-iy, 

i:,,sGaMp-q\\)G.P,-i) 

Since the local mode filtering is an iterative procedure which converges to the closest highest mode of the local 
histogram, one needs to propose an iterative scheme to solve the implicit equation 71 to get: 


” i:,esG<rMP-q\\)G>rXPq-Pp) ^ ^ 

for all p. 

So interestingly it has been observed that their is a good resemblance between bilateral filter equation with that of 
Equation 72. 

It has been observed that the bilateral filtering correlate to a gradient descent of a robust minimization problem. 
The problem definition: Input is a noisy image F. 

It is required to minimize the discrete equation: 


min - Ilf + Yj P(4 - 4))’ (74) 

^ peS qeN(p) 

where p and q are pixel positions and N{p) is the neighborhood of p, p is the weighting function. The second 
term of equation 73 is the regularization term which is sensitive to the intensity difference among neighboring pixels, 
influenced by the function p. This regularization term helps to find out the relationship with bilateral filter. To do that 
a new reweighed regularization term has been introduced to make the minimization problem take the form: 


mm GaX\\q-p\\)p{k-Ip)), 

peS qeN(p) 

To minimize equation 74 the following iteration is done: 


jt+i ^ 

” ” |A^(P)I 


G^,(||9-p||)p'(4-/p), 


qeN(p) 


By substituting p{s) = 1 - GcrX^) the equation becomes: 

A 


jt+i - f ^ 

" |A^(P)I 


E G^AIIq-pll)G<rP‘-il)(4-il), 


qeN(p) 


(75) 


(76) 


(77) 
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which is similar to the bilateral filter expression, leads to the weighted average as shown below. 

This form can be used for solving the minimization problem. Thus equation 76 and 77 clearly reveals the fact that 
bilateral filter is a robust filter. 

5.4. Robust Statistics and anisotropic diffusion 

Here we like to revisit the contribution of M J. black et al. ED who, with the very simple formulation showed the 
link between anisotropic diffusion and robust statistics. Equation 78 shows the continuous anisotropic diffusion filter 
form. 

^^^5gill=div[g{\\vl\\)vl] (79) 

Ot 

Similarly eqn. 79 shows the continuous form of the robust estimation problem which focuses on minimizing the 
equation as discussed in the previous section. 


nin r p(||V/||)JQ 

^ Jo 


(80) 


where Q is the image domain. Minimization of eqn. 79 is done via gradient descent using calculus of variations 
as shown: 


By substituting 


dl{x, y, t) 
dt 


= divlp (IIV/II) 


V/ 

m\^ 


(81) 


six) = 


p(x) 

X 


(82) 


the relation between robust estimation in eqn. 79 and anisotropic diffusion in eqn. 78 used for image reconstruction 
is obtained. 


5.5. Partial Differential Equation and Anisotropic Diffusion Algorithm 

In many imaging applications it has beeen observed that unknown data are related to observed data through a 
linear equation 

p = Rf + v (83) 

where v is the white Gaussian noise, is a linear operator defined on L^(7?^), considering the general problem of 
recovering unknown data / from noisy observed data p. Teboul et. al ISDl considered a minimisation function in 
order to estimate / from p defined as 

Jif)= f lf(x)-p(x)l^dx + A^ f <fi[lvf(x)l]dx (84) 

Jq Jq 

where Q is an open bound set in R^, |V/(x)| is the modulus of the gradient of f(x), A is a regularization parameter 
which balances the influence between the two terms of equation 81. Minima of (81), if they exist, formally verify the 
Euler Lagrange equation given by 


[f(x) - p(x)] - A div 


2^,,.^'[|V/(V)|] 


2|V/(x)| 


V /(x) = 0, xeEl 


dn 


bn = 0 


(85) 
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where n denotes the vector normal to the boundary JQ of Q, div is the divergence operator given by 


div{u) = 


2or3 ^ 
OUk 


The nonlinear equation (82) has a close resemblance to the anisotropic diffusion filter O. 
Equation 82 can be solved by a gradient-descent method given by: 


( 86 ) 


= [p(x) - fit, X)] A div— ———— V fit, v) = 0 


dt 


2 \Vfit, x)\ 


df 

-^IdQ = 0 
dn 


m 


where ns a scale parameter. The nonlinear PDE (84) corresponds to the anisotropic diffusion equation of Perona 
and Malik ||2l. 


6. Recent extensions/variations of bilateral filter 

There are several variants and extensions of the most extensively used bilateral filter. To name a few we like to 
discuss upon them in the next sequel. 


6.1. The Trilateral filter 

This trilateral filter as introduced by P.Choudhury et aLfTSl in the year 2003 is another variant and an extension 
of the bilateral filter for edge preserving smoothing operations. But it is meant for handling TV- dimensional signals 
in multimedia applications. Unlike their classical counterparts like anisotropic diffusion and bilateral filter, it is much 
robust for noise reduction and better outlier rejection in high-gradient regions. 

Applying the bilateral filter in regions of sharp change in gradient and high gradient areas degrade the smoothing 
capabilities of the bilateral filter as the rectangular filter extent encloses pixels spanning the peak of the ridge or valley 
as shown in Eigurej^ As a result the filter blends these intensities and results in a blunt feature. As shown in Eigure 
^the spatial extent of the filter can accommodate only a narrow portion of the input signal. These problems have 


been addressed by Choudhury et al. (13 by combining modified bilateral filters with a pyramid-based method to limit 
filter extent. 

Initially they have estimated the slopes of the intensity gradient using bilateral filter followed by tilting the filter extent 
by an angle 6 pivoting the center point at (x, fix)) to cover a wide area of usable domain. . Its tilting vector Ggix) 
should average together closely related neighborhood gradients ignoring nearby strongly dissimilar gradient outliers 
thereby assists in restoring the effectiveness of the spatial filter term. The tilting vector is given by: 


where 


G,(x) 


1 r 

k0(x) Jeo 


v/i„(x + Oc(Os(llv/i„(x + 0- v/i„(x)ll)d^ 


kff(x) 


%J oo 


ciOsi\\^Iinix + 0-^linix)\m. 


( 88 ) 


(89) 


c() denotes the domain of filter window, the range component. Instead of computing a range weight for 
neighboring fix + f) by measuring its closeness to the center point value fix), instead its closeness to a plane through 
fix) , which acts as a centerline for the filter window as shown in Eigure [3^ . It takes only one user set parameter, 
filters a signal in a single pass(thus non-iterative) as needed by most PDE based methods. As already mentioned it 
can handle A-dimensional signal and provides much better visual quality in applications like texture preserving in 
contrast reduction problem in digital photography and while denoising polygonal meshes. As shown in Eigure [33]\ 
and B. As mentioned by W.C.K Wong et al. ll4^ in their benchmark paper on medical image processing, that a narrow 
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Figure 32. (a,b and c) shows the unilateral, bilateral and trilateral filter range for one scan line of an image. Now given a (d) noisy signal (e) shows 
the bilateral filter smoothes the sharp corners(e.l) and averages the high gradient regions poorly(e.2) , (f) the trilateral filter smoothes corners(f.l) 
and high gradient regions(f.l) well, where 1,2 and 3 denotes ridge and valley like edges, high-gradient regions, and similar intensities in disjoint 
regions respectively. Figure reproduced from P.Choudhury El 


spatial window, say, 3 pixels in each dimension, should be used in order to avoid over-smoothing structures of sizes 
comparable to the image resolutions as sharp ridges and gutters which are commonly found in biomedical images, 
such as nested vessels in digital subtraction angiography (DSA) and 3D angiography, folded gray and white matters 
in brain MR imaging. Wide spatial window generally over-smooth sharp ridges, gutters and fine textures in an image 
E). So a trade-off is necessary to strike a balance between the size of the spatial window and the number of iterations 
needed to be performed in bilateral filtering. Their results have been shown in Figure 

T.Vaudrey et al. m in the year 2009 have shown that the trilateral filtering can further be speeded up several orders of 
magnitude (100-1000 times faster; from hours down to seconds) by the use of LUT(look Up Table) and by truncating 
the LUT to a user specified numerical accuracy. Referring to equation 82, cQ (the spatial weight) can be pre-calculated 
using a look-up-table (LUT) as t depends on the domain kernel standard deviation (say cri) and m (say the kernel size) 
and (^(offset from central pixel). N.B: size of the used filter kernel 2m-h 1 m is the half kernel size and is ^-dimensional. 
As this function is Gaussian, the LUT can be computed as: 


c(0 = exp(^) (90) 

where 0 < i < m. sQ cannot use a look up table since it is depends on the local neighborhood intensity. Numerical 
approximation needs to be done. The LUT based approach is a smart truncation of the kernel to a defined accuracy s, 
where 0 < £ < 1. If any value below s needs to be truncated then all the above it should be used given by, 

e < 

which leads to ||/|| < cri sl-2ln{s) = T 

where T is the threshold. This approach is not much beneficial to bilateral filtering but upon applying it to trilateral 
filtering, it reduces the number of computing equations drastically. This implementation is open to massive parallel 
processing (every pixel is independent within the iteration of trilateral filter) and opens up a huge opportunity for GPU 
(Graphics Processing Unit) based implementation. 

6.2. Joint and dual bilateral filter 

This filter commonly known as cross/joint bilateral filter is just a variant of the classical bilateral filter as introduced 
by Eisemann et al. ll44ll and Petschnigg et al. ll45]| . What it does: Given two input images I and E, the joint bilateral 
filter smooths /, while preserving the edge of the second image E. In other words the range component is computed 
using the image E given by: 
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Photographic Tone Trilateral Filter Original Model Noisy Model Smoothed Model 

Reproduction 


(A) 


(B) 


Figure 33. (A) In comparison with the bilateral filter, trilateral filtering limits blooming of sharp specular highlights and performs similar to the 
gradient attenuation method. (B) The single pass trilateral filter removes most visible corruptions caused by additive Gaussian noise in both vertex 
positions and normals. Figure reproduced from P.Choudhury im. 



(a) Original image (bjcloseup (c) Bilateral Filter (d) Trilateral filter 


Figure 34. MRI image (a) An original slice image, (b) magnified view of the slice, (c) After applying bilateral filter on the slice region, (d) After 
applying trilateral filter. Figure reproduced from Wong et.al. m. 
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Orig. (top) Detail Transfer (bottom) Flash No-Flash Detail Transfer with Denoismg 


Figure 35. The pic shows candlelit settings in flash and noflash environment. A particular area region of the image is magnifled as shown and the 
details of the image has been extracted with denoising as shown Figure reproduced from G. Petschnigg et al. EU 


JBFU, E]p=:^y G^XWp - ^\\)G<rXEp - Ep)I„ (92) 

'^P qeS 

with 

Wp = Y^ G^XWp - q\\)G,rXEp - Eq) (93) 

qeS 

and Wp defines the normalization weight, Gcr, and denotes the spatial and range kernel functions respectively. 
Figure [26| [27] and demonstrates the behaviour. 

Dual Bilateral filter: 

This is another variant of bilateral filter as introduced by Bennet et al. ll2^ . It takes two images as input I and J as 
input and both the input images are used to define the edges whereas he previous joint bilateral filter used only one of 
the hashed image to extract out the edge information. It is defined by: 

DBFmp =^y - q\\)Gcr,(Ip - Iq)Gaj(Jp - Jq)Iq, (94) 

qeS 

The merit of the above approach is that there remains very less probability to miss the edge information present in 
an image of various spectrum and in light or no-light condition. 

6.3. Nonlocal-means 

For each pixel p define a small, simple fixed size neighborhood as shown in Figj^ Secondly a vector Vp is 
defined which contains a set of neighbouring pixel values, i.e yp=[3 5 1 7 8 3...]. Similar pixels p and q\ have small 

vector distance(weight) given by: w{p,q\) = ||Vp - V^i|| and w{p,q2) = \\Vp - Vq 2 \\ . However dissimilar pixels 

such as p and q3 have large vector distances given by w{p, q3) = 11^^ - ^^311^- 
Therefore the non local means filter NL-means is given by 

NLMUlp \K - C (95) 

^P qeS 

where Wp = Y^qeS the normalizing factor. Vector distance to p sets weight for each pixel qi. 

Results shows a) is the noisy source image, b) is the figure obtained after gaussian filtering resulting in low 
noise and low detail, c) is the figure obtained after anisotropic diffusion filtering with some stairsteps as shown, d) 
is obtained after bilateral filtering which shows more clear picture than the previous ones but still some stairsteps 
remained which is prominent near the boundary edges, e) is obtained after applying the variant the NL-Means which 
exhibits sharpness, low noise and few artifacts. 
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Figure 36. Scheme of NL-means strategy. Similar pixel neighborhoods give a large weight, w(p,ql) and w(p,q2), while much different neighbor¬ 
hoods give a small weight w(p,q3). Figure reproduced from Antoni Buades et al. da. 




b 


Figure 37. a) Noisy input image, b) Gaussian smoothing, c) Anisotropic diffusion filter, d) Bilateral filter, e) NL-Means filtering. Figure reproduced 
from Antoni Buades et al. da 
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Figure 38. Illustrations of the bilateral filtering process (left) and the guided filtering process (right). Figure reproduced from Kaiming He et al. Eg. 


6.4. Guided filter 

Kaiming in their paper l|49l proposed a novel explicit image filter called guided filter. Derived from a local linear 
model, the guided filter computes the filtering output by considering the content of a guidance image, which can be the 
input image itself or another different image. The guided filter can be used as an edge-preserving smoothing operator 
like the popular bilateral filter lO, but it has better behaviors near edges. Some recent implementations of the bilateral 
filter fTl. lfTTll . llSOl have employed quantization techniques to accelerate the execution speed but at the cost of accuracy. 
The authors have shown unlike bilateral filter it does not suffer from the gradient reversal artifacts as shown in Fig 
It is defined as 

qi = akU + bk, 4iewk (96) 


where (uk, bk) are some linear coefficients assumed to be constant in Wk. The filtering output 
guided image. In order to determine the linear coefficients (uk, bk) constraints are needed from the 
The output q has been modeled from the input p by subtracting some unwanted components n like 


shown in fig 36 


= Pi 


is q and 1 is the 
filtering input p. 
noise/textures as 

(97) 


The cost function in the window Wk is minimized as 


E(ak, bk) = ^ (iakh + bk- pif + eap. 

iewk 

were 6 is a regularization parameter penalizing large ak. The solution of Uk and bk has been computed as 

^iewPiPi-pkPk 


bk = Pk- cikPk (100) 

where pk and o-f are mean and variance of I in Wk. Now Uk and bk being solved it can be used to compute the filter 
output q in equation 89. Besides the guided filter has also got some limitations like it may exhibit halos near some 
edges. Halos refer to the artifacts of unwanted smoothing of edges. It has been suggested by the authors to apply 
the guided filter in any situation when the bilateral filter works well. The guided filter is much faster and sometimes 
(though not always) works even better. 


(98) 


(99) 


7. Conclusions and future scope 

In this review article we have tried to present the benchmark edge preserving smoothing algorithms present keep¬ 
ing the main focus on anisotropic diffusion and bilateral filtering. We have discussed the various efficient and variant 
ways (with less computational complexity and time), to implement the edge preserving filters together with the inter¬ 
relationships with the structured numerical schemes with one another. We have discussed the various applications of 
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Figure 39. Detail enhancement. The parameters are r= 16, 6 = 0.1^ for the guided filter, and 16, cr^=0:1 for the bilateral filter. Figure reproduced 
from Kaiming He et al.(49l. BF= Bilateral filter, GF= Guided filter. 


the edge preserving filters keeping the focus on recent trends and their modifications and extensions in recent days as 
well, with an in-depth mathematical analysis. Our study on bilateral and anisotropic diffusion algorithm highlights 
several ways for future research. Moreover we have focused on hardware implementations of the most recent and 
efficient bilateral filter implementation and anisotropic diffusion filter with a view to make real time implementation 
making it a challenging task mainly in multidimensional data orientations used in computational photography. Future 
work will be to improve speed using parallel architecture (e.g GPU and using dedicated DSP blocks on reconfig- 
urable architectures) to meet the strict timing constraints keeping in mind the trade-off with energy efficiency. Our 
in-depth study on the edge preserving filters in this paper will help the future research community of the computer 
vision domain to attain a detailed evolutionary understanding of the edge preservers, from heat diffusion equation to 
the most resent trilateral filters, their advantages and disadvantages and their existing interrelationship, to evolve with 
new innovative ideas more efficiently. 
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